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Abstract 

When magnetic flux moves across layered or granular superconductor structures, 
the passage of vortices can take place along channels which develop finite voltage, 
while the rest of the material remains in the zero-voltage state. We present analytical 
studies of an example of such mixed dynamics: the row-switched (RS) states in un- 
derdamped two-dimensional Josephson arrays, driven by a uniform DC current under 
external magnetic field but neglecting self- fields. The governing equations are cast into 
a compact differential-algebraic system which describes the dynamics of an assembly 
of Josephson oscillators coupled through the mesh current. We carry out a formal per- 
turbation expansion, and obtain the DC and AC spatial distributions of the junction 
phases and induced circulating currents. We also estimate the interval of the driving 
current in which a given RS state is stable. All these analytical predictions compare 
well with our numerics. We then combine these results to deduce the parameter region 
(in the damping coefficient versus magnetic field plane) where RS states can exist. 

PACS numbers: 85.25.Cp, 46.10.+Z, 05.70.Ln, 47.54.+r. 

1 Introduction 

Two-dimensional (2D) arrays of Josephson junctions serve as "controlled laboratories" to 
investigate fundamental questions such as phase transitions U, vortex propagation and 
interaction 0, |3|, f|, [|, phase- and frequency-locking of coupled oscillators || 0, |8|, |9|, and 
spatio-temporal pattern formation and chaos [10, 11], among others |12|]. A standard circuit 



geometry is a rectangular array driven by a DC current uniformly injected from the bottom 
and extracted from the top in the presence of an applied field (Fig. p]). Their technological 
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Figure 1: 2D Josephson junction array consisting of N x = 7 rows and N y = 7 columns of 
square cells. The cell at is shown enlarged. Each junction is described by a gauge- 

invariant phase difference: <p x for the junctions on the horizontal edges, and (f) y for the 
vertical junctions. A uniform DC bias current I& c is injected into every node on the bottom 
edge and extracted from the top. The left and right sides are open boundaries. The mesh 
current ip denotes the deviation of the current distribution from a uniform current flow in 
the vertical direction. A uniform magnetic field /, in units of the flux quantum $o, is applied 
transversally to the plane of the array. 



promise as high-frequency oscillators [JT3], |TJ], [T5| depends critically on achieving tunable, 
highly nonlinear, coherent oscillations of the collection of junctions. However, to the chagrin 
of engineers (and to the curiosity of dynamicists) such coherent oscillations are not easy to 
obtain || |1^, |16| . Instead, the arrays frequently break up into incoherent substructures, and 
deliver output voltages with small AC amplitudes. 

A striking example of such dynamical states with spatial structure is provided by the 
row- switched (RS) states found in underdamped 2D arrays of square cells [0. As the bias 
current Jd c is ramped up, the DC current- volt age (i— V) characteristic of the array displays 
a succession of discontinuous jumps between ohmic branches of increasing resistance until, 
eventually, the normal resistive branch is reached. It was also noticed that the branches 
are equally spaced in voltage. These observations suggested a row-switching scenario, in 
which each of the jumps corresponds to all the junctions in an individual row suddenly 
switching from the superconducting to the normal state, thus increasing the voltage across 
the array by a fixed amount. In the RS states, the array then consists of superconducting 
and normal rows, coexisting to form striped patterns as in the four examples shown in Fig. [| 
In other words, the magnetic flux moves across the array along certain rows (channels) where 
a finite voltage develops, while the rest of the system remains in the zero-voltage state. This 
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row-switching picture was later explicitly confirmed by measuring voltages across individual 
rows [18[ and by direct imaging of the array [pOfl . Moreover, attempts were made [18| to 



determine the sequence of row-switching as the bias current was varied, and showed that the 
order and critical currents at which rows become normal are irregular and history-dependent. 

The row-switching phenomenon is robust and still appears when the underlying lattice 
structure is changed. In arrays of triangular cells it has been observed both experimen- 
tally plj and in simulations 



The explanation of row-switching could also be relevant 
to other systems. For instance, similar channeling of magnetic flux has been seen in contin- 
uous superconductors [p3 |. In addition, the hundreds of resistive steps which appear in the 
I-V characteristics of high-Tc superconductors |24j] have been taken as an indication of the 
layered weak-link structure in those materials. 

The experiments on 2D arrays of square cells have revealed that RS states appear only 
when the junctions are sufficiently underdamped [17, 19, 20. 25]. Otherwise, the I-V char- 



acteristics present an extended region of flux-flow leading to the ohmic branch of the entire 
array. In addition, RS states are only observed for sufficiently small applied magnetic fields. 
If the field is too large, there are no individual RS steps; rather, one giant step emerges |H| [T8| 
in the I-V. The origin of this giant step has been attributed to the interaction of self-fields 
with a coherent array oscillation in the form of a dynamical checkerboard pattern [26]. To 
understand such transitions between coherent and localized states, it is important to study 
the RS states and determine their current and phase distributions as well as the parameter 
regime for their appearance. 

Much of the previous theoretical work on 2D arrays has consisted of numerical simula- 
tions |27], ^8], |29], [3D], f|, |22|, |2|, |3_T], |32], [18], |25| which reproduce the measurements reasonably 
well. Several of them have discussed RS states briefly |29], |2| , while three others have treated 
in depth the row-switching phenomenon P2"| , [3^, 18]. Their main conclusion is the charac- 
terization of the two distinct types of rows found in the experiments: (1) switched rows ("S" 
rows), across which there is a finite DC voltage, and (2) quiescent rows ("Q" rows), across 
which there is no DC voltage drop. The simulations show that the junctions 4> y in the vertical 
branches of the 5* rows are in the normal resistive state (rapidly rotating) whereas those in 
the Q rows are nearly superconducting (stationary). Nevertheless, as we show, the junctions 
in the Q rows are still oscillating, which causes finite AC voltage drops and associated losses. 
This is why we hesitate to call the Q rows "superconducting" . 

The numerical investigations have also studied the sequence of row-switching as the bias 
current is varied. Even in the absence of temperature and disorder, the observed patterns 
and the order of their appearance are found [22, IS] to depend on factors such as the ini- 
tial condition, how currents are varied, the magnetic field (both externally applied and 
self- induced), among others. This is a clear indication that multiple attractors coexist for 
identical parameter values. (Indeed, Patterns 1,2, and 4 in Fig. |2| are found using the same 
set of parameters while Pattern 3 is obtained for a slightly smaller bias current.) When 
inhomogeneity is included, it becomes even harder to predict which row will switch next, 
except to conjecture that it will occur at the "weakest part" of the array [2~2~|. Phillips et 
al. have studied the vortex patterns in detail when inductances are included. When 
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Figure 2: Four snapshots of RS states in arrays of N x = 31 columns by N y = 7 rows. Two 
types of rows are observed: quiescent (Q) rows (in white) across which there are zero DC 
voltage drops, and switched (S) rows (shaded areas) across which there are finite DC voltage 
drops. Black dots denote topological vortices, defined in Sec. |[ They are (roughly) equally 
spaced in the S rows of the symmetric Patterns 1-3, but the spacing can change from row 
to row in asymmetric patterns such as Pattern 4. Correspondingly, their propagation speed 
(represented by the length of the arrows) may change from an S row to another within a 
pattern. In Patterns 1 and 2, the vortices move in phase, even when the S rows are separated 
by Q rows. These patterns are numerically generated using T = 0.2, / = 0.1 and I<j c = 0.6 
for Patterns 1,2 and 4, and Idc = 0.5 in the case of Pattern 3. Thus, Patterns 1,2 and 4 
correspond to coexisting dynamical attractors of the system. 
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self-fields are small, the S rows appear to be globally phase-locked even if they are far apart, 
separated by Q rows in between. This means that topological vortices in the S rows appear 
to propagate together, just as seen in Pattern 2 in Fig. |2|. However, for generic asymmetric 
patterns, such as Pattern 4, vortices do not move together. Stronger self-fields are also found 
p2| to break this phase coherence. 

Compared to the numerous experimental and numerical studies, analytical results are 
much scarcer for 2D arrays. As far as we are aware, previous authors have focused on 
the simplest solution, namely, when the whole array is on the normal branch of the I-V 
curve (Pattern 1 in Fig. This can be interpreted as the special RS state when all the 
rows have become normal; that is, the "completely row-switched" solution. These studies 
have concentrated on explaining the global phase-locking mechanism needed for oscillator 
applications. The complete RS state is found to be only neutrally stable under zero magnetic 
field [||, [33| (which implies that rows are decoupled), whereas a non-zero field induces inter- 
row locking. This inter- and intra-row coupling mechanisms have been studied by several 
methods: isolating two cells in the array |34], |7j, perturbation methods 0, and harmonic 
balance H . However, those results are not directly applicable to generic RS states since 
the complete RS state has no Q rows and it extends to any large bias current values for any 
damping. Instead, the (generic) RS states exhibit non-trivial structures and exist only in a 
certain parameter regime. 

In this paper, we study analytically the RS states and test our predictions against nu- 
merical integrations of the system. First, we cast the governing equations and the boundary 
conditions in the mesh formalism to ease the analytical procedure (Section [|). In this no- 
tation, the system can be viewed as an array of coupled oscillators in which the junction 
phases (the pendulum-like oscillators) are coupled via the mesh currents ip (the current 
distribution in the array). The coupling arises from the flux quantization condition. We ne- 
glect self-field effects in the equations, thus reducing the parameters of the system to three: 
the bias current Id c , the junction damping coefficient T, and the magnetic field /. In this 
way, many properties of the RS states can be explained without undue complications. We 
also discuss the notion of vorticity in these discrete arrays. 

We use primarily four examples (depicted in Fig. |2|) in order to illustrate and test our 
results. It is convenient to label each RS pattern by the set S of its switched rows. Therefore, 
Patterns 1 to 4 are labeled as: S = {1,2,3,4,5,6,7}, S = {2,4,6}, S = {4}, and S = 
{2,3,4,7}, respectively. We also define an S region to be a set of contiguous S rows. For 
example, Pattern 4 in Fig. |2| has two S regions, one with three rows 2-4 and another with a 
single row 7. Similarly, a Q region is a set of contiguous Q rows. 

A formal perturbation expansion in the high-frequency limit |35] is used to analyze the 
governing equations (Section ^). We assume that the RS states are time-periodic solutions 
in which some junctions whirl (i.e., the <f) y 's in the S rows are running oscillators), and all the 
other junctions librate (i.e., the <fi y1 s in the Q rows and all <p x are nearly stationary). Although 
the expansion is made systematic so that higher-order corrections could be obtained, we show 
that most of the features of the RS states can be accounted for by the leading order. (The 
only unresolved main feature is the phase- locking between S rows.) To the zeroth order, 
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we obtain two systems of algebraic equations: one for the DC, and another for the AC 
components of the phases and currents. The DC system is nonlinear (thus difficult to solve); 
however, we obtain bulk approximations which work well far from the edges. On the other 
hand, the AC components obey the linear discrete Poisson equation with forcing from the 
DC solution; therefore, they can be readily obtained once the DC solution is known. 

The bulk approximation determines analytically the DC and AC distributions of currents 
and phases for any given RS pattern. The first important result is that the DC current flows 
uniformly in the S rows, but circulating currents are induced in the Q regions. These strongly 
affect the spatial wave numbers of the S rows (also calculated analytically), thus explaining 
why the spacing and speed of propagation of the fluxoids in the S rows varies from pattern 
to pattern, and even from row to row within a pattern (Fig. |2|). In Section |] we test these 
findings numerically with good agreement. 

Another main conclusion from the leading order analysis is that the presence of S regions 
breaks the array into a collection of Q regions that are decoupled from each other, as far as 
the DC equations are concerned. The Q regions are, however, still weakly coupled through 
the AC component. Thus, for example, the existence of the switched row 4 in Pattern 3 
produces two 31 x 3 quasi-disjoint Q regions which only interact weakly. This picture proves 
useful because it reduces the problem of approximating the dynamical RS states of the array 
to obtaining the static states of its (smaller) constitutive Q regions. 

Indeed, this physical picture has further implications for the stability of the RS patterns 
(Section |5|). As explained above, each RS state is only observed in an interval of the bias 
current, which depends on the magnetic field and damping. We show that the upper current 
limit of this interval is well predicted by the depinning current of the largest Q region. This 
means that the RS state ceases to exist when the flux penetrates any of the Q regions which, 
in absence of irregularities, is usually the largest one in the array. For example, Pattern 3 
cannot hold beyond the current value at which a static state of a 31 x 3 array depins. Along 
the same lines, we also argue that the largest upper current of any RS state will correspond 
to RS patterns whose largest Q region is a single row, such as Pattern 2. However, the 
depinning currents are independent of the junction damping, whereas the RS states are 
found only in underdamped arrays. This indicates that the lower current limit also plays an 
important physical role. A crude approximation for this lower limit is the retrapping current 
of a single junction which does depend on the damping. Combining these criteria we are 
then able to calculate the region of the parameter plane of the magnetic field / versus the 
damping parameter Y where RS states are not possible. Throughout Section |] we present 
additional numerical evidence which support these criteria. 
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2 Formulation 



There are two equivalent ways of formulating the governing equations of the system: the 
node and mesh formulations. The node formulation is easier for simple geometries but 
it becomes cumbersome and impractical for two-dimensional arrays when inductances are 
included. Thus, we follow the previous literature []27|, |29], [|, [31], |32], [l8j , and derive a compact 



description of the arrays in the mesh formulation. In particular, we follow closely Phillips 



et al. [plj |32|1 and Trias [18], with a few changes. Although this formalism was originally 



developed to ease numerical simulations, it is well suited for analytical work. 
2.1 Governing equations 

Our description of the array shown in Fig. [I] assumes several simplifications. First, we 
neglect thermal fluctuations (i.e., zero temperature), and we consider all junctions identical 
(i.e., no disorder). Second, we describe our basic circuit unit, a single Josephson junction, 
by the resistively and capacitively shunted junction (RCSJ) model. In this standard model, 
a junction driven by a current I b is represented by an equivalent circuit of three channels in 
parallel with a capacitance C, a resistance R, and a tunnel junction with the critical current 
I c . As a result, its state variable (the gauge-invariant phase difference across the junction) 
is governed by 

M[4>\ =0 + r0 + sin0 = / 6 (1) 

where the nonlinear operator M returns the total current through the device. In ([J) the 
current is normalized by I c , whereas time is expressed in units of the inverse of the plasma 
frequency lu' 1 = (%C/2nI c ) 1 / 2 . In addition, T = f3~ 1 / 2 = (%/2-Kl c R 2 Cfl 2 is the damp- 
ing, with f3 c the McCumber parameter. Also, $o is the quantum of magnetic flux. The 
instantaneous voltage across the junction is given by the Josephson voltage-phase relation: 

V(t) = T0, (2) 

where the voltage is normalized by I C R. Thus, a single junction is analogous to a damped 
driven mechanical pendulum, and its voltage corresponds to the rotation frequency of the 
pendulum |35], |36], ^7 . 



When several junctions are interconnected to form a network, like the one depicted in 
Figure [I], the current distribution must fulfill Kirchhoff 's current law. This results in coupling 
among the junctions. It is convenient to decompose each branch current into an external 
and a deviation current: 

I b = hxt + Idev (3) 

The external current J ex t is chosen such that it satisfies current conservation at all nodes, 
including external sources and sinks. In general, it can be spatially non-uniform or time- 
dependent. However, in our case, the steady bias current I^c is injected (extracted) at the 
nodes along the bottom (top) edges. Therefore, our choice [[J^] for J ext is the stationary 
uniform vertical flow, in which I ext = I<i c on every vertical branch of the circuit (for all t), 
and /ext = on every horizontal bond. 
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The deviation I dev from the external flow must be divergence-free since current sources 
and sinks have already been incorporated. Therefore, there is a stream function (or mesh 
current) ip at each cell whose discrete curl determines Jdev in the x and y directions: 

r dev (hj) = Hhj)-Hhj-i) (4) 

i y dcv (hj) = -m,j)-^-i,j)). (5) 

(In the rest of this paper we will not write time-dependences explicitly when they are obvious, 
such as here.) 

In order to ensure that these relations hold also at the edges of the array, we define 
artificial boundary cells which have either the horizontal index % = or N x + 1, or the 
vertical index j = or N y + 1. This yields the boundary conditions of the problem: 

ij)(i,j) = ifi = 0,N x + l or if j = 0,W + 1. (6) 

This condition is equivalent to "grounding" the value of ip outside the array. 
Combining we obtain the first two sets of governing equations 

M[<i> x (iJ)] = ^(ij)-^(ij-l) (7) 
MW&j)] = h c -(^(i,j)-^(i-i,j)) (8) 

where Af was defined in ([j]) . 

The other source of intrinsic coupling between the junctions is due to a macroscopic 
quantum constraint: the flux quantization condition around a cell. Each corner of a cell is 
a superconducting island described by a well-defined phase. Calculating the phase change 
around cell yields the third and final set of equations of the system 

m% + l,j) - - (4> x (i,j + 1) - r(hj)) + 2^%^ = 27rn(z,j) (9) 

for i — 1, . . . , N x and j = 1, . . . , A^^, where is the total magnetic field penetrating 

the cell. The winding numbers n(i,j) are a set of integers that arise because the island 
phases are only defined up to multiples of 2tt. The initial condition sets n(i,j) which remain 
constant as long as the array is kept superconducting. However, without loss of generality, 
all n(i,j) can be set to zero. Suppose they are not zero; then we can redefine the junction 
phases as 

o" (/../) - o'(i.j) 
y (z,j)-27r£n(A;,j) - <jP(i,j) (10) 

k=l 

such that (H) is unchanged except, now, alln(i,j) = 0. Crucially, both currents and voltages 
are invariant under this redefinition of the phases since adding integer multiples of 2tt to 
<j) y changes neither sin0 y nor cf) y . This means that the dynamics and measurements remain 
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identical for any combination of integers n(i,j), and we do not need to be concerned with 
their initial values. Similarly, if the magnetic field were controllable independently on each 
cell, adding an integer number of flux quanta $o into any cell would not change the measured 
I—V characteristics, at least within this model. This is simply the array analogue of the 2- 
junction SQUID, whose dependence on the penetrating field is also ^-periodic. Because of 
this periodicity in the magnetic field, the topological vortex must be defined differently in 
2D arrays and in continuous superconductors, as we will discuss at the end of this section. 
The total magnetic field in (|||) can be decomposed into two parts: 

$(l,j)= $ ert + $tad(«,j)- (11) 

The first term is produced by the external field applied perpendicularly to the plane of 
the array, which we assume to be constant and uniform. It is usually parametrized as a 
dimensionless frustration / normalized to the flux quantum: 

/ = <W$o, (12) 

such that, in terms of /, the period of the external field is unity. The second term, the 
induced field, can be written generally as the sum of all the contributions from the branch 
currents 

$ind(U)=EE^4 b (13) 

n k 

where k runs through all the branches of the circuit while n corresponds to the four edges 
of cell The branch inductances L b nk are purely geometric constants determined from 



the circuit. PL 31 



2.2 Matrix-vector notation 



The above equations can be cast into a compact matrix- vector notation ||31|| . For a N x x N y 
array, all branch variables (e.g., currents I b , voltages V, and phases 0) can be written 
as vectors of dimension equal to the number of branches, i.e. (N x + l)N y + N x (N y + 1). 
Thus, for instance, the vector consists of all the phases <p x and <p y . On the other hand, 
variables defined at cells (e.g., the mesh current ip and the induced flux ^Va) form vectors of 
dimension N x N y . These two groups of vectors are connected via a branch-to-cell connectivity 



matrix |39[ M which takes a directed sum as we loop around a cell: 



Mcf>(i,j) = + 

- + (14) 

More mathematically, this operator takes the discrete curl of <fi around every cell 
Conversely, the discrete curl of the cell variables is obtained by applying the transpose cell- 
to-branch matrix M T . 

Using this notation, the total flux (|TTD can be written as 

$ = $ Q / + ML h I b . (15) 
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where L b is the branch inductance matrix, and / is a constant vector. 
Moreover, (£D,fl5|) can now be written simply as 



/ dev = M^, (16) 

and become 

N[<t>\ = hxt + M T ip (17) 

where Af is operated component-wise and the vector J ext has components that are zero on 
the horizontal edges and Idc on the vertical edges, as defined by our choice of uniform vertical 
flow. 

Finally, we can use flI3|) and (|TJ) to recast the flux quantization condition (£5|) as 



M0 + 2nf + ^-(L m ij + ML b I cxt ) = (18) 

where components of L b are normalized to flop, p is the lattice constant, Aj_ = $o/2vr/ c /i p 
is the dimensionless penetration depth, the mesh inductance matrix is defined as 

L m = MlJ'M 1 . (19) 

and we have set n(i,j) = 0. 

To summarize, the governing equations (|lj]) and (fL8"D form a closed differential-algebraic 
system for <p and ip, with parameters /, T, I^ c , Aj_, and the coefficient matrix L b . This form 
of the system is compact and intuitive. It can be seen as a coupled-oscillator system in which 



the "oscillators" <p are driven by the "coupling field" ip in (17). In return, the oscillators 
collectively feed back onto the field in (|18"D. This picture suggests the following integration 
steps [^2[ [L8||: first, given <p at some time t, solve fll8|) for L m tp\ then, invert the matrix 
L m , together with the boundary conditions ([]), to determine the field tp. This gives us the 
"drive" on the right hand side of fllTD, which is used to update the oscillators <p in time. 



We conclude the general formulation by pointing out that the equations flTT|), fli"8] ) possess 



two simple symmetries p9| . If we find a solution ((p(i,j,t),ip(i,j,t)) for / and 1^, then 
(—<fi(i,j,t),—ip(i,j,t)) is a solution of the system for — / and —Idc, the other parameters 
being the same. Similarly, (—<p(—i,—j,t),ip(—i,—j,t)) is also a solution for / and —Idc 
(since M is changed to —M due to the reversal of the spatial coordinates). Therefore, we 
only have to study the quadrant / > and Jd c > 0. Together with the unit periodicity in 
/, the frustration can be further restricted to0</<l/2 without loss of generality. In the 
rest of this article, expressions such as "large /" and "small /" are used within this interval. 



2.3 No-inductance approximation 

Computing the full equations QT7D,(|T8|) quickly becomes a heavy task as the system size 
increases. In previous numerical studies, these computational limitations have been circum- 
vented either by using acceleration schemes 0, when the inductance effects are of interest 
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per se, or by "truncating" the matrix L m (i.e. neglecting some of its components). Three 



truncations ||, |31], [32], [TjJ are often used: no-, self-, and nearest-neighbor-inductances. Self- 
inductance neglects the inter-cell magnetic coupling by keeping only the diagonal components 
of L m (which then becomes trivially invertible). Nearest-neighbor- inductance includes, in 
addition, magnetic coupling between neighboring cells. An important remark is that not only 
the mesh inductance L m but also the vector ML b I ext must be provided in order to complete 
the system, and the choice of J ext may affect the results when L b is truncated j|, [18|]. (In 
contrast, the choice of J e xt is unrestricted if the full inductance matrix is used.) Truncating 
the system in a physically consistent manner is a subtle issue, and, for simplicity, we shall 
assume no inductance in this article. 

In contrast to what one might guess from its name, the no-inductance approximation 
does not eliminate the inter-cell coupling. Counterintuitively, it leads to an even longer-range 
coupling than the self- and nearest-neighbor- truncations. The no-inductance approximation 
sets L b = 0, thus L m = 0. The flux quantization condition flTE| ) is then just 

M0 + 2nf = 0. (20) 

The same equation can also be obtained in the limit Aj_ = oo for any L b , which allows the no- 
inductance limit to be approached from the inductive system continuously. It is important 
to note that the condition (|20| ) is now a constraint on <f), which must be satisfied at all times. 
The discrepancies between M0 and —2nf cannot be filled by locally adjusting the induced 
field, which was possible when the inductive terms were present. This leads to a global 
coupling of the junctions over the whole domain. To see the coupling mechanism provided 
by (pUt), operate the loop-sum M on (|T~7|) . From the left hand side we obtain 

MN((j))(i,j) = M4> + TMcj) + M[sin0] 

but the first two terms vanish since fl2"0"| ) must hold identically. From the right hand side, we 
obtain 

M(J ext + M Tr ip) = MJ ext + MM' c = M/^ - 
where we have introduced the discrete Laplacian 

Aijj{i,j) = + -1) 

+ iP(i + l,j)+i>(i-l,j))-4i>(i,j). (21) 

For the stationary uniform flow J ext , the term MI ext = 0. Thus, we arrive at a discrete 
Poisson equation 

= — M[sin </>] (22) 

in which the distribution of the mesh current is dependent on all the junctions in the array. 

Equations ( |TTD and ( |2"2"| ) constitute the governing equations for the no-inductance case, 
and can be integrated in the same manner as before. Provided that the initial condition 
satisfies the constraint ( PU| ) and its time derivative M0 = 0, (^) is satisfied for all t. An 
immediate advantage of the no-inductance approximation is that the sweep of the parameter 
space is greatly simplified since the number of parameters has been reduced to three: /, T, 
and /de- 
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2.4 Vorticities in 2D arrays 

Before closing this section, we consider now the concept of vorticity in these arrays. Analo- 
gously to incompressible planar fluid flows, we can define a vorticity by taking the curl (by 
applying M) of the "velocity" field which, in our case, corresponds |39| to the branch current 
I b . This current vorticity 

n = MI b = M(J ext + M T tjj) = MJ ext - A^j (23) 

measures how currents whirl, and can take any real values. (For the uniform vertical external 
flow the contribution MI ext vanishes.) 

In contrast, the notion of a topological vortex (or charge) is commonly used in Josephson 
arrays in analogy to continous superconductors. In Type-II superconductors, the vortices 
would correspond to the integer winding numbers n(i,j) in the flux quantization condi- 
tion (g). But, as we showed above, the n(i,j) are irrelevant in the arrays. Therefore, an 
alternative, less physical definition for the topological vorticity has been used: 0, (|, |32], |1^ 



C = i_ M (0-0). (24) 

Here, <fr denotes restriction of the components of the phase vector <fr within [— ir, n). The value 
of ( at each cell takes only integer values (typically or ±1) and jumps discontinuously as 
the system evolves in time. In effect, this definition detects when one of the four junctions in 
a cell rotates and crosses (f> = n (mod 27r), since M<fi changes discontinuously by 2tc at that 
instant. This is the 2D analog of marking the location of a ID kink at the point where = 7r 
(mod 2ir) regardless of whether the kink really has a localized structure or not. This particle- 
like picture is frequently useful but, by neglecting spatial distributions, there is a potential 
loss of information about the true dynamical state of the system. On the other hand, the 
current vorticity Q would reveal more accurately how localized vortices are. However, our 
simulations show that the topological vorticity ( moves together with a peak of the current 



vorticity f2 (Sec. [4.2|) . Thus, we use both definitions interchangeably at our convenience. 
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3 Analysis 



In this section we present a perturbative analysis of the governing equations (|T7|) , (p0|) . Al- 
though the analysis is made systematic so that it is possible to proceed to higher orders, we 
show that most of the fundamental features of the row-switched states can be explained by 
the leading order of the expansion. 

Before writing down the RS solutions, it is useful to think of the array with uncoupled 
junctions. This limiting case corresponds to imposing tp — in ([l?]), thus reducing the 
array to a collection of uncoupled pendula, independently responding to a constant drive. 
The junctions on the horizontal branches (0), whose drive is zero, converge asymptotically to 
4> x * = due to the damping. On the other hand, the uncoupled vertical junctions (|8]), driven 
by Jdc, have a different dynamical behavior. For small T, they can converge asymptotically 
to one of two distinct stable states P6[ : the superconducting (static) state, which exists 
only when I& c < 1, in which the drive is balanced by the sinusoidal nonlinearity, (i.e. <p y * = 
arcsin/dc); or the ohmic (whirling) solution, where the first time derivative balances the 
drive, and increases at a nearly uniform rate u = idc/T (i.e. the pendulum "whirls"). The 
two attractors may coexist for the same drive, and hysteresis may occur. 

When the junctions are coupled, the simple dynamics of the independent junctions is 
altered, and complex spatio-temporal solutions, which do not have an analogue in the un- 
coupled array, may emerge. Nevertheless, in the case of the RS solutions the two states of 
the driven single junction mentioned above (static and whirling) are still valuable "building 
blocks" for the analysis of the whole system. Specifically, the RS states are characterized 
by alternating regions in which the vertical junctions are either stationary (Q regions) or 
whirling (S regions). There are, however, some significant differences with the uncoupled 
case. For instance, the time-averaged current distribution in the coupled array deviates from 
the uniform flow. Hence, the phases of the stationary junctions can have other values than 
or arcsin/dc- In addition, the rotations of the vertical whirling junctions induce AC os- 
cillations on the stationary junctions, and phase-locking among the whirling junctions. Our 
analysis in this section is capable of explaining most of these effects. 

We note that our analysis is restricted to solutions with no (static) vortices trapped in 
any of the Q regions. The "no-vortex" state is expected to be most relevant to determine 
the parameter regime in which RS states appear. Similarly, the vertical junctions in the S 
regions are assumed to be whirling at almost constant frequency neglecting more nonlinear 
rotation, which is certainly possible, as we discuss in Sec. 

3.1 Perturbative analysis 

In previous perturbative analyses of junctions and arrays, it has been customary to treat 
large parameter [^, [|(J. However, partially RS states can exist only when J dc is 



sufficiently small, as we will show below. Therefore, we use the rotation frequency of the 
pendulum to = Jdc/T as the large parameter in our perturbation. That is, we will consider 
the high-frequency limit [35] u ^> 1, which can be satisfied for a finite Jdc if the damping T 
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is small enough. 

Hence, we assume that the variables in the RS states can be expanded in powers of u^ 1 . 
The phases of the horizontal junctions are then approximated by 

oo 

a, (i,j,t)=W(i,j) + Ew-^(i,j,r). (25) 

p=2 

while the mesh current is given by 

oo 

ip(hj,t) = ipo(i,j) +Vo(i,J,r) + J2 u ~ P ^p(hji T ) ( 26 ) 

p=\ 

where we have introduced the normalized time 

T = Ujt = (Jde/r) t. 

The notation (•) expresses time-independent (DC) quantities, while (■) are for the time- 
dependent (AC) parts whose time average is zero. Note that the correction of 0(cj _1 ) 
in (^) turns out to be zero, so we neglect that term from start. The form for the vertical 
junctions must be different in the switched and the quiescent rows. In the switched rows, all 
the junctions are whirling and their phases grow, to the lowest order, constantly in time: 

oo 

^(i,j,t) = T + $ i (t,j) + J2^ y P (hJ,T), jeS. (27) 

p=2 

Meanwhile, in the quiescent rows the junctions are librating and, thus, the leading order is 
stationary: 

oo 

^(U,f)=^(i,^ + £0^(1,^7-), J'eQ. (28) 

p=2 

We impose ipo and the higher order terms to be periodic in time. In general, the period has 
to be modulated and, thus, expanded in u^ 1 (strained coordinate). However, since in the 
following we will focus on the leading order system, we set the period to be exactly 2tt in r 
for simplicity. 

The perturbative calculation proceeds in the usual way by substituting (|2l)|)-(|2"B|) into 
(|17]) , (^(]) ; Taylor-expanding the sine in ([!]); and equating terms of the same order in to. In 
principle, this procedure can be carried out to higher orders if secular terms are eliminated 
by satisfying solvability conditions when they arise. 

Balancing the leading order terms, we obtain two sets of equations since the time- 
independent (DC) and time-dependent (AC) terms must cancel separately. First, the DC 
terms yield the following equations for both types of rows: 

sin0g(i,j) = if) Q (i,j) -ipo(i,j ~ 1) (29) 
+ <^(z + l,j)-<^(z,j) (30) 
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and one more equation which depends on the type of row (switched or quiescent): 



l dc -sm j) 



: iJ>o{hj) - 4>o{i - j e S 



(31) 
(32) 



These equations constitute the full DC system. 

Similarly, from the AC terms we obtain for both rows 

<!>f(ij,i~) = ip (hj,r) - if) (i,j - l,r) 
4>l(i,j + l,T) - (P x 2 (i,j,r) = 



(33) 



(34) 



where " denotes differentiation twice with respect to r. Moreover, for each type of row we 
obtain a different equation: 



and 



02 = -sm(T + (/) y (i,j)) 

- ^o(hj,r) + ip (i - l,j,r), jeS 

<ii\hj,T) = -i/j (i,j,r) +ip (i - 1, j,r), j G Q, 



which completes the full AC system. 

These systems of equations are to be solved with boundary conditions 

= tpo = 

at the boundary cells. 

A simple but important observation can be made at this point. Using 
i = and N x + 1 (i.e., at the right and left edges), it follows that 

^o(i,j) = ° Vz, if jGS\ 



ID and 



(35) 
(36) 

(37) 
at 
(38) 



Therefore, the leading order DC mesh current vanishes in a switched row [43], just as it does 
in the top and bottom boundary cells at j = and N y + 1. In other words, the switched row 
is equivalent to having another boundary row, which splits the array into two. Thus, to the 
leading order, a partially row-switched array with many switched rows can be described as 
a collection of disjoint quiescent regions, coupled only weakly through the AC component. 
This useful picture is exploited later. 

The solution of the leading order systems is otherwise non-trivial since the DC equations 
(p9f)-(|32]) constitute a nonlinear algebraic system, and the DC solution is in turn needed 
to solve the AC system (p3|)-(p6|). Thus, in general, they have to be solved numerically - 
although we show below that useful approximations can be obtained under certain assump- 
tions. 
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Once the leading order solutions are found, the calculation could be carried out to higher 
orders. The next order correction leads to a particularly simple set of equations: 

4>f{hj) + T<tf(i,j) = MiJ) ~ MiJ - 1) (39) 
4>l"(i,j) +T4> y 2 '(i,j) = -Mhfi+M* - U) (40) 
<%(i,j + 1) - <%{i,j) = $(« + 1,3) ~ <f>l(hj) (41) 

for all r and regardless of whether the row j is switched or quiescent. Terms from the 
sinusoidal nonlinearity do not come into play at this order, but further expansions would 
certainly involve more complications. 

It is important to note, however, that the salient features of the solutions observed in 
the numerics can be explained from the leading order equations. Therefore, we restrict our 
analysis to the DC and AC equations in the following sections. On the other hand, we will 
also point out a remaining problem which is likely to be resolved only by considering the 
higher order corrections. 



3.2 Analysis of the DC equations 

The DC equations (J29f)-(|3"2f) constitute a nonlinear algebraic system which must be solved 
numerically in general. However, to gain insight into the system, we will now obtain approx- 
imate solutions to the system when there is a large asymmetry between its two dimensions. 
We will then come back to the full system and discuss its solutions. 



3.2.1 Large aspect-ratio approximation 

Consider the case when all quiescent regions in the array are longer horizontally than verti- 
cally. This happens, of course, when the array itself satisfies N x ^> N y . More importantly, 
arrays whose dimensions do not fulfill this condition are also broken into smaller, laterally- 
long, almost disjoint quiescent regions after several row-switching events. Thus, this "large 
aspect-ratio" approximation is important to characterize the RS states which appear in the 
course of the row-switching process. Remember we also assume that none of the Q regions 
contains static vortices, which could be trapped for large N y and /, and for small I& c . In Sec- 
tion |5| we will give an estimate of the values of / and N y for which we expect this assumption 
to be valid. 

In such a situation we expect a nearly "uniform" solution in the bulk of the array with 
some edge corrections near the right and left boundaries. Hence, far from the boundaries, 
we assume the vertical junctions in the quiescent rows to become independent of i, 

<t>o(hj) = 4>l(j) forjeQ. 



On the other hand, we assume a whirling solution for the switched rows in which waves 
with well-defined wavenumbers k(j) propagate: 

Mh3)~-k(j)i + S(j) for j E S. (42) 
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Note that the wavenumber k(j) and the phase constant S(j) may differ from one switched 
row to another. The other DC variables 0g, and V'o are also assumed to be i-independent. 
Thus, the DC equations reduce to 

sin0!(j)=TO)-^(j-i), (43) 



r (j + J L) - 0g(j) = W - k(j) for j G 5, (44) 
<f%U + 1) - 4> x oU) = 2tt/ forj eg, (45) 

W) = forj65, (46) 

0o(j) = arcsin J dc for j G Q. (47) 

This simplified set of equations is still nonlinear but solvable. We begin by analyzing all 
the quiescent regions in the array (if any), delimited by switched regions or by the physical 
boundaries. Consider a quiescent region spanning from row ji to j 2 (> ji). Such a region 
contains n = j 2 — j\ + 2 rows of horizontal junctions including the top and bottom borders, 
and n — 1 quiescent rows of vertical junctions. We emphasize that these vertical phases are 



all given by (|47|), thus, I& c < 1 is necessary for the existence of partially RS states, where Q 



rows are present. From fl43|) the horizontal phases must satisfy a telescope sum 

J2 + 1 

J2 sin0g(j) = MJ2 + 1) - MJi " 1) = 0, (48) 

where we have used the fact that both rows j 2 + 1 and ji — 1 must be either switched or in 
the boundary cells, and thus ip = from (f46|) or (|37|) . Now, ( |45|) can be solved with ( [481) to 



obtain: |42| 



0§(j + h - 1) = 2tt/ - —J (49) 

with j = 1, . . . ,n. This gives the time-averaged phases of the horizontal junctions in the 
bulk of the Q region. Then, from (^), the mesh current in the same region can be computed 

as 

MJ + Ji ~ 1) = E sin m) = — T^r sin [nf(j - n)] (50) 

for j = 1, ... , to — 1. This procedure allows us to solve for each Q region in the array 
independently. 

The remaining variables are easy to find. Recall that ipo vanishes everywhere in the S 
region. The rest of the horizontal junctions 0g lie either between two S rows, or between an 
S row and a boundary cell. In either case, it follows from (E3[) that 



0o 0) = 0; inside a S region. 

Finally, the wavenumbers k(j) for the switched rows (j G S) can be calculated from 
j). One notices that k(j) can change from a row to another, depending on the adjacent 
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horizontal junctions 0g. On the other hand, if there is an S region with three or more rows, 
the wavenumber k(j) = 2irf for all the rows except for the two rows at the top and bottom 
borders of the region; this is because 4>o = inside the region. In this sense, k = 2ixj is the 
"natural" wavenumber for S rows. 

This concludes the solution of the simplified equations (]4"3|)-(|4"7D . We now exemplify this 
procedure with four RS states of an array with N y = 7 rows, as depicted in Fig. |2|. In 
Section f| we will compare the predictions with our numerics. 

Pattern 1: S = {1,...,7}. 
This is the totally row-switched state in which there is no Q region. Thus, the horizontal 
phases are 

0g = (0,0,0,0,0,0,0,0). 

(The j-th component of the vector is 4>o(j)- Note this j runs through 1, . . . , 8 for the 7- row 
array.) In addition, ipo(j) = 0, and k(j) = 2i\f = k for all rows j = 1, . . . , 7. 

Pattern 2: S = {2,4,6}, (and so, Q = {1,3,5,7}). 
In this symmetric pattern there are four Q regions, each consisting of only one row, and 
three one-row S regions. By solving each Q region independently, we find 

0g = 7r/(-l, 1,-1, 1,-1, 1,-1,1). 

Then, for the S rows j = 2,4,6, we have ipo(j) = — sin(7r/) and k(j) = Anf. That is, the 
three S rows have an identical wavenumber but different from the natural ko- 

Pattern 3: S = {4}. 

In this case the two symmetric Q regions, rows 1-3 and 5-7, are separated by the central 
S row. We obtain: 

0§ = vr/(-3,-l,l,3,-3,-l,l,3). 
The wavenumber of the S row is fc(4) = Sirf. 

Pattern 4: S = {2,3,4,7}. 
In this highly asymmetric switching pattern there are two Q regions. We obtain: 

0g = 7r/(-l, 1,0, 0,-2, 0,2,0). 

The S rows have the following wavenumbers: k(2) = 3irf, k(3) = 2irf, and k(4) = k(7) = 
An f. Note that the rows 2-4 are contiguous but all have different wavenumbers. The row 
3 is surrounded by other S rows, hence has the natural wavenumber. Meanwhile, the rows 
2 and 4, which are contiguous to Q regions have different wavenumbers. 

A similar bulk approximation can be obtained for the other limit of the aspect-ratio. We 
present this small aspect-ratio case in Appendix 

One might wonder what has happened to the phase constants 5(j) of the switched 
rows ([42"D. Indeed, the 8(j) have disappeared in the simplified system (f43])-((47|), making 
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them arbitrary. However, simulations show that the switched rows are weakly coupled, so 
that the 8's drift to some particular values (if / ^ 0). This phase-locking has been noticed 
in the completely switched state and left unexplained |8|, [RJ. As we show in the numerics 
of Section |], it is also a feature of the partially RS states. The indeterminacy of 5 in our 
analysis is not merely due to the assumption of the whirling solution (f42|). Rather, it is 
already inherent in the DC equations (|29|)-(|32"|) for which the addition of a constant to all 
the </>o (i) within any switched row leaves the system unchanged. Since the drift occurs in a 
much slower time scale than the basic oscillation frequencies [33|, we conjecture that the S(j) 



could be determined from solvability (or secularity) conditions that might arise from higher 



orders of the expansion. That was the case in one-dimensional series arrays [[HJ where a 
similar slow phase drift and eventual locking was explained in that manner. However, it is 
beyond the scope of this paper to develop a similar calculation for the 2D array, and, in the 
following, we will use the values of S(j) obtained from the numerical simulations. 

3.2.2 Solving the full DC equations 

We now consider how to solve the full DC system beyond the bulk approximation — a problem 
which requires, in general, numerical solution. It is important, however, to note that the 
decoupling of the equations introduced by the switched rows is still present so that the 
problem reduces to calculating static solutions of smaller arrays. 

The important point to recall is given by (^): the mesh current is still zero in all S 
rows. This breaks the array into disjoint Q regions, as far as the leading order DC part 
is concerned. Mathematically, this means that equations (p9|),(|30D, and (j32|) are closed 
within each Q region, and can be solved independently. This system is identical to the 
superconducting (static) equations for an isolated 2D array of the same size as the Q region. 
When this sub-problem of finding the static solutions for the independent Q regions is solved, 
the remaining unknowns, 4> V Q in the 5* regions, can be determined from (|30|) . This two-step 
procedure is completely analogous to the one used in the large aspect-ratio approximation, 
except that 4> v in the Q rows now depends on i, and, thus, 4>o(i,j) in the S rows cannot have 
the form given in (0). 

How do we obtain the static configurations? Since a Q region can take any size in the 
j-direction (up to N y ), we need, in short, a general calculation scheme of static states for 
an arbitrary rectangular array. An analytical formula is not known even for the no-vortex 
solutions (one of the many possible superconducting states) we are primarily concerned with. 
Thus, they must be found numerically JGJ. A rare exception is the ladder array, of size 



N x x 1, for which an accurate analytical approximation has been obtained j|7]]. It shows that 
the full static solution differs from the bulk approximation in the existence of skin layers near 
the left and right edges. Crucially, resolving the phases in the skin layers is central to the 
existence and stability of the static solution. The ladder case is special but important since 



it is the most persistent in the parameter space among Q regions of a given width |E| [16 
In Section [5] we will connect the stability of the RS patterns with the stability of the static 
states. 
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3.3 Analysis of the AC equations 



We now study the AC system (p^)-(|36|). We only need to note that this is a lin ear system 
which is forced by the sinusoidal drive sin(r + 4>q). Therefore, if the DC solution is known, 
the AC system is simple to analyze. 

Assuming that the homogeneous part simply decays, the solution locks to the forcing and 
the time-dependence can be factored out as: 



" 01 " 




' A ' 






B 






C 



(i, j)exp(r v /r T) + c.c. 



(51) 



where "c.c." denotes complex conjugate. Then, the spatially-dependent complex amplitudes 
must satisfy 

-A(iJ)=C(iJ)-C(iJ-l) (52) 
-B{i,j) = -C{i,j) + C{i - 1, j) + f(ij) (53) 
A(i,j + 1) - A(i,j) = B(i + 1, j) - B(i,j) (54) 

with 

f{iJ) = [ ^xp j)V=T) if j e S 

I o if j G Q- 

Eliminating A and B from the equations, we obtain a discrete Poisson equation for C: 

AC = -n (56) 

with the source term 

»(i,j)=f(hj)-f(i + hj) (57) 
and, from (|37D, boundary conditions 

C = in the boundary cells. (58) 

In the rectangular domain this problem can be solved via the double discrete Fourier-sine 
series 
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2(Ny + 1) 



(61) 



2(N X + 1), 

Finally, A and B are determined from (|52"D and (0). This completes the analysis of the 
leading-order equations. 
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4 Numerics 



4.1 Finding RS states in simulations 

To test the validity of the analysis developed in the previous section we now compare its 
predictions with numerical results. The full governing equations (|I7|) and (pOf) , together 
with the boundary conditions (||), are integrated using the standard fourth and fifth order 
Runge-Kutta integrator with adaptive time step. Ours is an elementary non-optimized 
version of the previous mesh-formulated code [TBI] which enables us to switch between 



no-inductance and simple inductance models. The results presented here are all obtained 
neglecting inductances. 

As stated above, we also neglect the effects of temperature and disorder. Since most 
of the analysis have assumed the large aspect-ratio approximation, we study an array with 
N x = 31 and N y = 7, with small damping T = 0.2 and a moderate external field / = 0.1. 
We use as initial conditions the predicted large aspect-ratio DC approximations 0q ,s/ (and 
the corresponding first time derivatives). They are expected to be close enough to the true 
RS states to facilitate convergence, but we leave the AC part to be adjusted by the system. 
We choose a value for idc between and 1, and monitor whether the ensuing dynamical state 
is indeed the attempted RS pattern. 

The system, of course, does not always converge to the row-switched state we have 
targeted; the chosen initial condition may be out of the basin of attraction of the target 
state, or the state may not exist, or it may be unstable for the chosen parameters. The 
outcome from using "wrong" parameters is, as far as we have tested, as follows: If Idc is too 
large, then vortices start to enter in some of the rows we have initially set quiescent; if idc 
is too small, then the rows we have set switched cannot maintain the whirling motion, and 
exhibit retrapping, become quasi-periodic, or show highly nonlinear oscillations. In those 
cases, we then adjust idc until we find the clean periodic RS solutions which we aimed at. 

Not only must we tune Jd c , but the damping parameter T must be small enough in order 
to find clean RS states. If T is too large, it is difficult to find any partially RS states at all. 
For intermediate values, such as T = 0.4, some RS patterns are observed, but some others 
cannot be found. For the underdamped case T = 0.2 studied here, it becomes easy to find 
an appropriate range of Idc i n which the system converges to the expected RS pattern. This 



dependence on the damping is in qualitative agreement with experimental findings [17 . It is 
also consistent with our assumption of the high-frequency limit since a smaller T for a given 
Idc corresponds to a larger uj. 

Generally, patterns with large quiescent regions are more difficult to obtain; for example, 
the RS state S — {1} (with one Q region of 6 rows) has a smaller interval of suitable Jd c 
than the symmetric Pattern 3 S = {4} (with two Q regions of 3 rows), even though both 
states have only one S row. These observations and the above parameter dependences will 
be discussed in more detail in Sec. [|. 

Before presenting detailed comparisons between numerics and analysis for Patterns 1-4, 
we first illustrate convergence in Fig. |3|. There we show the time evolution of two variables 
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Figure 3: Time evolution in the central (switched) row of Pattern 3: (a) horizontal phase 
4> x (16,4,t) and (b) mesh current ^(16,4, t). After a short transient, the solution converges 
onto a periodic attractor with DC values and AC amplitudes (shown together as the bands 
delimited by the dotted lines) well predicted by the analytical formulas in Sec. §. 
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Figure 4: Time evolution of the vorticities in the middle of the array of Pattern 3. The 
solid curve depicts the current vorticity Q(16, 4, t), while the topological vorticity ((16,4, t) 
switches discontinuously between (no vortex) and 1 (one vortex in the cell). This dis- 
continuous "tagging" of the position of the vortex is clarified by the dotted curve, which 
corresponds to cos0 y (16, 4, t). Inspection of that magnitude indicates that everytime it be- 
comes —1 (i.e., the phase is equal to ir) one topological vortex enters the cell (and ( is 
increased by one). 

in the array for Pattern 3, using I^ c = 0.5. Since the initial condition (taken as the bulk 
approximation) is not a solution of the full system, there is a short transient (t < 50) until the 
system settles onto a periodic attractor. Recall that only row 4 is switched in this pattern. 
Figure |3](a) shows the phase X (16,4) of a horizontal junction adjacent to the switched row 
and in the middle of the row, where the large aspect-ratio (bulk) approximation is expected 
to be valid. The approximated average value is 2>nf ~ 0.94, as predicted in Sec. |[ Similarly, 
the mesh current in the central cell ^(16,4), shown in Fig. |3|(b), is ip = on average with 
some oscillations, as expected in any switched row. Not only the average values but the AC 
amplitudes are also well estimated from the AC leading order equations, as demonstrated in 
Fig. |. 

4.2 Vortex motion 

We illustrate now the two vorticities defined in Section [| In Fig. |4] we show the current 
vorticity Q and the topological vorticity ( both at the central cell (15,4) of Pattern 3 after 
convergence. They display similar periodic behavior though Q is continuous whereas ( 
switches discontinuously between and 1. ( becomes unity when a charge enters the cell, 
which occurs in this case when <f) y (lG, 4), the left junction, crosses n (modulo 2n). Therefore, 
( becomes unity when cos0 y (16,4) = —1, as shown in the figure. Similarly, when the right 
junction 0^(17,4) (not shown) turns and crosses tt, the charge ( is reset to zero. 

As a complement to the time evolution of the vorticities in one cell, we now show snap- 
shots of their spatial distributions for all Patterns 1-4 in Fig. |5|. Each cell is shaded according 
to the value of the current vorticity j): dark regions indicate positive large Q, while bright 
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Figure 5: Snapshots of Patterns 1-4 showing spatial distributions of the current vorticity 
fl as density plots. Dark regions correspond to large positive Q. Compare them with Fig. |2] 
where the same spatial patterns are shown in terms of the topological charge (. We observe 
that the topological vortices are generally located on peaks of fl, and propagate locked to 
the underlying wave. 
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parts correspond to negative f2. The same snapshots, but showing the topological charges £, 
are given in Fig. |[ Even though Q represents the spatial structure more clearly, we observe 
that a charge in a cell corresponds to a peak of Q, and that the charges propagate through 
the array on top of the underlying wave. 

Thus, we can use ( to visualize the wavelength and the propagation speed in each row. 
In all Patterns 1-4, the charges move across the array at a nearly constant speed, as seen 
in the space-time plots of ( in Fig. [| They propagate only through the S rows, and are 
apparently in-phase in all rows for Patterns 1 and 2. However, in Pattern 4 the S rows 
are not in-phase, and the propagation velocities vary from row to row. Thus, the simplistic 
picture that vortices carry all the flux and move with the same speed in all the S rows 
within a pattern leads to estimated speeds in disagreement with our simulations. This 
further proves that the underlying assumption that the topological vortices are particle- 
like objects which concentrate the flux is not accurate. Instead, the RS solutions are not 
localized states and the flux is spatially distributed, as suggested by previous work |TD[ 
and demonstrated in our analysis. Therefore, in these states, the topological vortices merely 
mark where the rotating junctions cross tt (mod 27r) (see Sec. |2.4| ), and they travel at the 
phase velocities of the underlying (non-localized) waves. Our analysis correctly estimates 
the spatial wavenumbers (thus, the propagation speeds) as shown below. 



4.3 Spatial structures after convergence 

We now present a quantitatively comparison of the analysis of Sec. ^|to numerical simulations. 
The analytical predictions correspond to the large aspect-ratio (bulk) approximation both 
for the DC and the AC components. For the numerics, we simulate a 31 x 7 array, and the 
system is allowed to converge to periodic solutions for Patterns 1-4 using T = 0.2, / = 0.1 
and Jdc = 0.6 (except for Pattern 3, in which I^c = 0.5 had to be used). 

We first check the predicted spatial wavenumbers k(j) in the S rows just discussed above. 
In Fig. [?| we show a "snapshot" of the (j) y in the S rows (2,3,4,7) of the non-trivial Pattern 4. 
To ease the display and comparison of the numerical results, we have juxtaposed the rows 
one after the other. Within each row, the spatial dependence is clearly linear, thus justifying 
the whirling mode assumption ([421) . The predicted wavenumbers &2 = 37r/, k% = 2irf, and 
&4 = &7 = Airf (dashed lines) are almost indistinguishable from the numerics (solid lines) 
except for small deviations close to the edges. 

Recall that in our analysis of the DC equations the inter-row phase differences 6(j) are 
predicted to be arbitrary in (|42|). Hence, only the slope of the spatial dependence is known 
and the dashed lines are adjusted to match at the center of each row. Conversely, this is 
a way to determine the 6(j) from the numerical simulations. For the four Patterns, we 
obtain: 

Pattern 1: 5(1) = 6(7) = 0, 6(2) = 6(6) = 0.05, 
6(3) = 6(A) = 6(5) = 0.1. 

Pattern 2: 6(2) = 6(A) = 6(6) = 0. 



25 




Figure 6: Space-time plots of the propagation of topological vortices for Patterns 1-4. The 
vertical (space) axis is the cell index: the cell is indexed one-dimensionally by i + 

N y (j — 1) by juxtaposing row after row. Within each of the symmetric Patterns (1-3), the 
vortices in the switched rows have the same wavelength and are in-phase. However, in the 
asymmetric Pattern 4, the spatial wavenumbers differ from row to row. 
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Figure 7: A snapshot of the vertical phases <p y in the switched rows S = {2, 3, 4, 7} of Pattern 
4 in a N x = 31 by iV y = 7 array. Each solid line connects the numerical phases of the 32 
junctions in each switched row. The dashed lines (almost overlapping with the solid ones) 
are the analytical approximation in Section ^| which predict the observed spatial wavenumber 
very well. The horizontal axis denotes the "vertical edge index," which numbers the vertical 
junctions consecutively as i + (N x + l)(j — 1) for i = 1, ... , N x + 1 and j = 1, . . . , N y . This 
enables us to display the 2D array in a single axis by juxtaposing one row after the other. 
As a guide to the eye, vertical dotted lines are added to separate the rows. 



Pattern 3: 5(4) = 0. 

Pattern 4: 5(2) = -1.8, 5(3) = -4.7, 5(4) = 0.2, 
5(7) = 0. 

Note that in each case one S(j) is set to zero and taken as the reference, which is equivalent 
to choosing the origin of t. In the following, we will use these numerical values of 5 when 
needed (most importantly, for the analytical values of the AC components). 

Next, we compare the predicted DC values with the numerical mean values after conver- 
gence in Figures |8| and |9|. As we showed in Fig. |](a), each horizontal junction <fr x librates 
around some DC value after convergence. These average values are plotted (solid lines) and 
compared to the large aspect-ratio approximation (dotted lines) in Fig. The prediction 
is uniform within each row because edge effects were neglected — consequently, it works well 
everywhere except close to the right and left ends. In Fig. ||| we show the DC values of the 
vertical junctions <p y (i,j). Note that the phases of whirling junctions are unbounded, thus 
not shown. In particular, Pattern 1 is left out from this figure since all <p y are switched. 
Again, in the bulk of the array the approximation (<j) y = arcsin/dc) holds, but near the edges 
there is a significant deviation. 



In a similar manner, we show in Figs. [LQ and 11 the AC amplitudes of the horizontal and 



vertical junctions, respectively; that is, the |A(i,j)| and \B(i,j) \ calculated in Section Ol As 



seen in both figures, symmetric Patterns 1-3 have rather constant amplitudes throughout 
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Figure 8: The DC values of <p x for Patterns 1-4, showing the spatial distribution of the 
average horizontal phases. The horizontal axis is the "horizontal edge index," defined as 
% + N x (j — 1) for <f) x (i,j). There are N y + 1 = 8 horizontal edges so that j runs from 1 to 8. 
For each j, the N x = 31 phases in the same row are connected. The dotted lines are from 
the large aspect-ratio approximation which accurately estimates the numerical results in the 
bulk of the array. The DC values are predicted to be multiples of irf. The approximation 
neglects the effects of the left and right edges, and, thus, inevitably misses the skin layers at 
both lateral boundaries. Vertical dashed lines mark the separation between fs. 
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Figure 9: The DC values of (j) y for Patterns 2-4. Switched junctions are unbounded, thus 
not shown. Therefore, Pattern 1 is absent in this figure since all the junctions are switched. 
The horizontal axis is the vertical edge index, defined in Figure |7|. For each quiescent row 
its N x + 1 = 32 phases are connected. The large aspect-ratio approximation, dashed lines 
at <j) y = arcsin/dc; is a good approximation in the bulk, but misses the lateral skin layers. 
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Figure 10: Dimensionless AC voltage amplitudes (f> x for Patterns 1-4, plotted versus the hor- 
izontal edge index. The large aspect-ratio approximation is shown as dashed curves. There 
are some quantitative discrepancies, but the approximation captures the spatial distribution 
of the amplitudes. 
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Figure 11: Dimensionless AC voltage amplitudes cjfl for Patterns 1-4, plotted against the 
vertical edge index. Again, the large aspect-ratio approximation, shown as dashed curves, 
can describe the spatial distribution fairly well. 
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each row, except near the left and right edges. The asymmetric Pattern 4 shows spatial 
fluctuations. Our estimates, shown as dotted lines, reproduce the spatial structure fairly 
well. It is quite remarkable that our approximation roughly captures the behavior at the 
right and left boundaries, when we have used the bulk approximation <J)q (together with the 
numerical to solve the AC system. 

Since the mesh current is determined from the phase configurations <j) x ' y , it also compares 
well with the large aspect-ratio approximation. Thus, we do not display the quantitative 
comparison of ijj, and instead present more descriptive 2D contour plots of the numerical ip 
on the 31 x 7 array geometry. The contour curves of the DC component of ip are shown in 
Fig. [12|. If the 2D array were continuous, the induced currents would flow along these curves 
on average. Since the array is discrete, the flow is restricted to the branches, but the level 
curves still describe roughly the way the currents circulate. Furthermore, the DC values in 
the S rows are nearly zero, as expected. In the Q regions, currents circulate in the clock-wise 
direction (ip Q < 0) on average. This would induce a magnetic field through the Q regions in 
the opposite direction to the external field /. Although it is interesting to ask whether the 
induced field cancels the external one to produce a Meissner-like region, that question only 
makes sense when all inductances are included. Note also that all contour plots are almost 
left-right symmetric, but show a slight asymmetry. This is presumably due to the presence 
of edges and the preferred direction (x) of propagation of the waves across the S rows. Such 
details are not captured by the bulk approximation and the full solution of the DC equations 
becomes necessary. 

Finally, the amplitudes of the AC oscillations of ip are shown in another set of contour 
plots in Fig. |T^. The observed nodal structures (typical in linear forced systems in a bounded 
domain) show the spatial distribution of the modes locked to the driving DC solution. The 
magnitude of these AC amplitudes is comparable to the DC values in Fig. |T^, even though 
the oscillating components of the phases <fi are much smaller, of O(uo~ 2 ), than their DC 
values. This is consistent with our analysis which assumes that the mesh current has DC 
and AC components both of 0(1). 
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Figure 12: Level curves of the DC mesh current ip for Patterns 1-4, indicating how the 
induced circulating currents flow. The total current flow is the superposition of the induced 
flow and the injected uniform current flow. Contour levels at —0.1, —0.3, . . . , —1.1 are drawn 
on the 2D grid of the N x = 31 by N y = 7 array. Pattern 1 shows little deviation from the 
uniform current flow on average, thus ip = and no curves appear. In the other patterns, the 
DC values of ip in the switched rows are zero, while the values are negative in the quiescent 
rows. Therefore, currents circulate in the clockwise direction in each quiescent region "along" 
the level curves shown. Strictly, the currents are restricted to the grid, but the level curves 
provide an intuitive description of the flow. Note that the boundary condition ip = is 
imposed at a half cell outside of the array borders; this explains why some of the contour 
curves intersect the array edges. 
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Figure 13: Level curves of the AC amplitudes of the mesh current ip for Patterns 1 (top) 
to 4 (bottom), on the 2D grid of 31 x 7 cells. Contour levels at 0.1, 0.2, . . . , 1.5 are shown. 
The magnitudes are generally large in the switched rows, but even quiescent rows have 
some oscillations and, thus, are not purely superconducting. Our leading-order analysis in 
Sec. |37| predicts that these AC oscillations obey the discrete Poisson equation with forcing 
originating from the DC components. The figure shows nodal structures typical in solutions 
to such a problem. 
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5 Parameter region for RS states 



In this section we determine where in the parameter space we expect RS solutions. This is 
a difficult task, partly because the parameter space is large. Even after neglecting induced 
fields (i.e. Aj_ = oo) we are left with three parameters: /, T, and idc- In addition, there can 
be multiple attractors coexisting for a given parameter set. Recall, for example, how in the 
previous section Patterns 1, 2 and 4 were obtained using an identical parameter set, and 
Pattern 3 also used a similar Idc value. A thorough determination of the parameter regime 
would then require a rigorous study of the bifurcations of the branches of all the attractors — 
an exploration which exceeds the scope of this article, and is perhaps too detailed to justify 
the necessary effort. Here, we take a more heuristic approach, and make several assumptions 
to estimate the current interval [J min , J max ] in which a given RS state is an attractor, as a 
function of / and T. We base our assumptions on the results of previous sections, and we 
demonstrate their validity by additional calculations in the following. 

5.1 Upper current limit 

We first estimate the upper current J max at which a given RS state ceases to be an attractor. 
Our first assumption states that this upper limit is reached when vortices enter any of the 
Q regions from the edge. The entrance of flux might produce further switching of rows 
(resulting in another RS state where the original Q region has been subdivided), or a more 
complicated state where the flux remains static or moves through the original Q region in a 
highly nonlinear motion. In either case, the original RS state is no longer maintained. As 
discussed in Sec. [5], each Q region is decoupled up to the DC leading order and is equivalent 
to an isolated superconducting array of the same dimensions. If, as we assume, no vortex 
has been trapped beforehand in the Q regions, arrays with more rows depin at smaller values 
of Jdc, as can be shown numerically pEEfl . Therefore, our second assumption is that, as Jd c is 
raised in an RS state, a vortex first enters the largest of the remaining Q regions, causing 
further break-up of the array. 

Thus, once the depinning current values for no- vortex superconducting state of any num- 
ber of rows is known, these two assumptions enable us to estimate the upper / max limit for 
any given pattern. For example, Pattern 3 (S = {4}) has two Q regions of the same size (3 
rows). We expect then that this state is not sustainable beyond the depinning current of a 
31x3 array. At zero temperature and without disorder, the likely scenario is that flux enters 
the center row of each of the two regions, so that a new RS state, Pattern 2 (S = {2, 4, 6}), 
ensues. This state has now four Q regions, each consisting of one row. The upper I dc value 
for this state should coincide with the depinning current of the 31 x 1 "ladder" array. Be- 
yond this value all rows switch and Pattern 1 is obtained. We have indeed observed such 
a sequence of row-switching events when we gradually increase Jdc from zero, using a clean 
initial condition: = = everywhere. Similarly, the largest Q region in Pattern 4 has 
2 rows. Therefore, Jdc should coincide with the depinning current of a superconducting no- 
vortex 31x2 array. In Table [I] we summarize the excellent quantitative agreement between 
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Table 1: Stability intervals [J mm , J max ] (two middle columns) for eight RS patterns (two of 
them identical) in the 31x7 array using / = 0.05 and T = 0.2. The set S denotes the switched 
row numbers, and Patterns from Fig. |2] are labeled. The intervals are calculated numerically 
by gradually changing Jd c and following the corresponding branch of the RS state until 
instabilities appear. For example, Pattern 2 is found in the interval [0.335, 0.945], for this 
set of parameters (/ = 0.05, T = 0.2). The upper limit 7 max can be predicted accurately by 
the depinning current 7 dep of the largest Q region of each pattern (with dimensions N x x N y , 
shown in parentheses). The lower limit I mm is harder to estimate, but the retrapping current 
F et of a single junction serves as a rough estimate: for V = 0.2, the value is F et = 0.252, 
which is smaller than the observed I mm = 0.305-0.335. The first four rows show Patterns 
1-4 from Figure |2|. The next four patterns all have a single S row, but its location is different. 
Among these four, Pattern 3 has the widest stability interval because its largest Q region 
(31 x 3) has the smallest number of rows. 

the numerically observed / max values of several RS patterns, and the depinning currents of 
superconducting arrays with the same dimensions as their largest Q region. 

We have also tested our assumptions with four additional patterns, all with only one 
switched row: S = {4} (the symmetric Pattern 3), S = {3}, S = {2}, and 5* = {1} (the 
most asymmetric pattern). This illustrates the dependence of the upper J max not on the 
number of switched rows, as above, but on their location. For given / and T, J max becomes 
smaller as the switched row is shifted from the middle of the array to the bottom because 
the largest Q region increases its size from 31 x 3 to 31 x 6. Excellent agreement is again 
obtained between our criterion and the numerical observations (Table [I]). 

We now make the third assumption that enables us to estimate / max analytically in some 
cases. We propose that, as the drive increases, just before the entrance of a vortex into a Q 
region, a junction barely holds itself at a critical angle 

crit = ±7r/2. (62) 
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When it is forced to turn beyond that value, depinning takes place, just as it would if the 
junction were uncoupled. Recall that the single uncoupled junction under an increasing 
drive becomes unstable through a saddle-node bifurcation at Jd c = 1, with <fi = 717 '2 as the 
bifurcation angle. Although the criterion for global depinning is different in a coupled array, 
this simple heuristic criterion has been used to predict the depinning current in ladder arrays 



with remarkable accuracy f47H 



Take, for instance, an array at zero temperature with small N y in a ground state with 
no pre-trapped vortices. Then, the first junction to cross <p crit = ±7r/2 is, for / > and 
Jdc > 0, the vertical junction who sits in the center row at the left edge. Thus, the flux 
would penetrate the array through that junction and destroy the RS state. This is readily 
deduced from the circulating current shown in Sec. |] which reinforces the drive near the left 
boundary. Such a current is due to the presence of the left and right boundaries, which our 
large aspect-ratio approximation neglected. A full analysis of the skin layers would be needed 
for a general analytical prediction, but there are two tractable limiting cases of interest. 

The first case is a "small aspect-ratio" superconducting region, i.e. with many more rows 
than columns (N y ^> N x ) . As discussed in Appendix [A], a bulk approximation can then 
be used, which approximates accurately the phases near the left and right edges — because, 
in this case, the skin layers are located near the top and bottom boundaries. Should such 
a region be present in a RS state as a Q region, it would be very easily broken even with a 
small value of Jdc- In fact, J max can be quickly estimated as the drive for which the central 
leftmost vertical junction crosses the critical angle fl6"2"|). As shown in fl6"8|), this junction 
has indeed the largest angle. The depinning value for such a small aspect-ratio region is 
estimated to be 

J sar- 2{NX + 1) jl+ |- W 

From (|63|), the region remains stationary when < ic$^ and / < 1/2N X . If / > 1/2N X , 
a vortex enters the N y ^> N x region for any Jdc > 0. We have tested these conclusions 
numerically with good agreement. Moreover, note that other physical arguments |37| predict 



that the edge barrier for the penetration of flux in this limit would be roughly given by 
f c ~ \/txN x . The condition (|63"D results from the instability of a static state, and it does not 
depend on T, the damping coefficient. 

The second case is the "ladder array", with N x columns and a single row |[48| . Its 
superconducting states, including states with trapped vortices, and their bifurcations have 
been studied comprehensively [|47[| . One of the results of that work is the curve of the 
depinning current as a function of /, shown as a solid line in Fig. |14|(a) and (b). This 
monotonically decreasing curve is again independent of T, and becomes insensitive to N x , as 
soon as N x is greater than about 5. For / up to about 0.46, the depinning is caused by the 
disappearance of the no-vortex solution. This part of the curve is well approximated by the 
solution of the following implicit nonlinear equation which comes from imposing ( |62"D to 
the leftmost junction: 

arcsin(l - 1^ ) + ^ arccos(/ L m A - ) = tt/ (64) 
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Figure 14: (a) Stability region for Pattern 2 with T = 0.4. For / = 0.05, 0.1, 0.2, we sweep 
/dc to determine numerically the stability interval [J mm , J max ]. denoted by the vertical arrows 
with endpoints at J mm (o) and J max (•). The solid line is our estimate for J max (/), given 
by the depinning current J dcp of a ladder array. The dashed line is an estimate for 7 mm , 
given by the retrapping current J re ^ of the single junction at V = 0.4. Therefore, the shaded 
section is the estimated region of the /dc - / plane where Pattern 2 exists, for T = 0.4. Note 
that the region does not extend beyond a critical / = /rs(0.4). (b) Same as (a) but for 
T = 0.2. Although the upper estimate J dep is unchanged, the lower estimate 7 re ^ decreases 
with T. Consequently, Pattern 2 is expected to be observed in a larger parameter region for 
smaller T, as shown by the five intervals (arrows) obtained numerically. The region does not 
extend for / larger than / RS (0.2). (c) Phase diagram for the existence of RS states in the 
f-T parameter plane. The curve / = /rs(F) separates the regions in which no RS states 
may or may not appear. For r > T* m 1.2, J re ^ = 1, thus, no RS states are expected for any 
/. This diagram explains the previous (qualitative) observation that RS states occur only 
when junctions are underdamped and the applied magnetic field is small. 
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with r = a + \/a 2 — 1 and a = 1 + y 1 — /lad 2 / cosnf. After a crossover at / ~ 0.46, the 
static checkerboard pattern becomes more robust, and this formula ceases to be valid. We 
will not discuss RS states in this high range of /. If our assumptions are correct, this critical 
curve should predict the J max of Pattern 2. In addition to the single comparison presented 
in Table [I] for this pattern, we test it for T = 0.2 and 0.4 and several values of / in Fig. |14[ 
As shown there, the numerical J max values of Pattern 2 from simulations are predicted very 
accurately by the analysis of the ladder depinning point. 

Up to now, we have assumed that the magnetic flux penetrates the Q regions from the 
left edge of the array. However, the flux can also enter the array from the top or bottom 
boundaries of a Q region in certain situations. Consider a Q region with a large aspect-ratio 
and no trapped vortices, but when the number of rows, say N y , is large. In this case, the 
bulk approximation obtained in Sec. [3.2.1| can still be used. From ((49"|) the maximum angle 
for the horizontal junctions is <p x = N y irf attained at the top and bottom edges of the region. 
It is clear that this value becomes larger than the critical angle (|62|) when / > l/2N y . Thus, 
for a fixed N y while / is increased, the flux would enter the Q region roughly above that 
value of the frustration. The entrance of flux in this manner puts a limit on the applicability 
of our analysis. The assumed no-vortex Q region is expected to exist only when the number 
of rows is smaller than about 1/2/. Thus, our analysis does not apply for the initial stages 
of the row-switching cascade in large arrays, when there are still Q regions with many rows. 
However, even in such arrays, later steps of the cascade (when the Q regions have been 
subdivided) can be described by assuming no- vortex Q regions. In addition, our preliminary 
simulations indicate that keeping a vortex trapped in a Q region becomes more difficult both 
in the presence of Idc (which tends to expel the fluxoids from the Q regions) , and of self-fields 
(which tend to shield the Q regions from the entrance of vortices). 

5.2 Lower current limit 

The parameter regime for the existence of RS states is also affected by the damping parameter 
T, as shown in experiments and simulations where RS states appear only when the junctions 
are underdamped. However, all the critical currents calculated until now are independent of 
T. We claim that the explanation of the T dependence of the RS states requires an estimate 
of the lower limit 7 min . Unlike the upper limit, in which the superconducting solution of a Q 
region ceases to exist, our numerical observations suggest that the lower limit is caused by 
an instability mechanism in an S region. As the bias current is decreased from the values in 
which a clean periodic RS state is observed, S regions start to have trouble in maintaining 
fast whirling oscillations. Typically, the system begins to show amplitude modulations in a 
slow time-scale, becomes highly nonlinear, or gets retrapped altogether. 

The variety of possible scenarios makes an accurate estimate much harder than for / max . 
In order to make progress, we have to rely on a rather rough estimate, based on the dynamics 
of a single junction. Recall that the vertical junctions in an S row are in the resistive 
(whirling) state. For a single underdamped junction, its inertia is enough to maintain a 
whirling solution until very close to the retrapping current 7 ret , when it jumps back to the 
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stationary state. Only near that value does a strong nonlinearity come into play. Ignoring 
the inter-junction coupling, we use this current as our estimate for the lower limit 7 min of an 
RS state. Because of collective effects, the state may not be immediately retrapped into a 
stationary state, or, on the contrary, be retrapped earlier. However, we expect that, as Ia c is 
lowered toward the J ret value, some nonlinear effects start to become apparent, so that the 
simple periodic RS state is altered. 

The estimate of J ret is standard For the underdamped case, (i.e., r < T* « 1.2), 



the retrapping is produced through a homoclinic bifurcation at I vct < 1, and the I—V of the 
single junction is hysteretic. For all T > T*, F et = 1, and there is no hysteresis. In general, 
7 ret is calculated numerically, but an asymptotic expression, J ret ~ AY/n can be used as 
r -> 0. 

From the definition of J ret , our estimate for J min is thus independent of /, N x , N y , and 
the particular RS pattern, but depends on the damping T. The estimates for T = 0.2 and 

0. 4 are shown as dashed lines in Fig. 0(a) and (b), respectively. The comparison with the 
numerical values of J min (the point when the RS states lose their whirling character) is not 
so good, as expected. However, our estimate seems to serve as a reasonable first guess. 

5.3 f—T parameter region for RS states 

In the usual experimental setup, the I—V characteristic of an array is measured by sweeping 
the DC current under a constant applied magnetic field at a fixed temperature (which controls 
the penetration depth Aj_ and damping T). For some combinations of the experimental 
variables (magnetic field and temperature) and, thus, of the underlying parameters /, T, 
and Aj_, the I—V shows RS steps. For others, it does not. We will now summarize the 
preceding sections and combine their results to estimate the (T,f) parameter region, in the 
limit A_i_ = oo, in which RS states appear. 

First, in Sec. |5.1| we showed two limiting cases in which J max can be obtained analytically, 

1. e. when N y = 1 (ladder) and when N y ^> N x . Numerical simulations |22], [47], |46|] show 
that J max changes monotonically between these two limits, as N y is varied. (This result is 
also expected from physical grounds: for fixed N x , the magnetic flux penetrates the array 
more easily as N y is increased.) An obvious consequence of this is that the ladder array has 
the largest parameter domain for the no-vortex superconducting state. Therefore, recalling 
our link between depinning and row-switching, RS states whose Q regions are all ladders, 
e.g. Pattern 2, are thought to be the most stable in the same sense. In other words, when 
an isolated ladder of length N x cannot maintain superconductivity, the 2D array of size 
N x x N y cannot show row-switched behavior. Consequently, the solid curve in Fig. |T4](a) 
not only gives the upper limit J max for Pattern 2, but also establishes the critical Jd c , for 
each /, above which no RS states can be observed. 



Second, we concluded in Section [5]2] that the I mm of all RS states with damping T can 
be estimated by the retrapping current of a single junction with the same T. These are the 
dashed straight lines in Fig. |H|(a),(b). 



Hence, the RS states can only exist in the region contained between these upper and 
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lower limits, shown as the shaded area in Fig. 11(a). Those limits intersect at a value /Rs(r) 
beyond which no RS state is possible. (See Fig. [I4|(a)-(b) for the procedure.) For junctions 
of moderate to large damping (r > T* « 1.2), the dashed line is above the curve, meaning 
that RS states are impossible for any /. On the other hand, for highly underdamped arrays 
(r < 0.2), the line always remains below the curve; hence, RS states are possible for any 
/ (although the region of / near 1/2 would need more careful consideration). Between 
these two extremes of damping, the line intersects the curve at the critical value /rs(F), 
which constitutes a phase boundary in the f-T plane. In other words, the parameter plane 
is divided into two regions (RS and no-RS) by the curve /rs(F) in Fig. 0(c). This is in 
qualitative agreement with previous observations, and awaits more systematic experimental 
testing. 
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6 Summary and open problems 



In this article we have used a weakly-nonlinear perturbative analysis to study the row- 
switching phenomenon and to approximate the RS solutions. For the bulk of the array we 
have obtained analytical expressions for the phase and current variables. In addition, we 
have estimated the parameter regime for their existence. For this, the consideration of the 
lateral edges has played an important role. The predicted spatial current distributions and 
the parameter regime could serve as a guide for more systematic experimental studies. In 
the rest of this section we briefly state open problems and possible future directions. 

The leading-order solutions show good agreement with the numerics, but leave one phase 
per row undetermined. This is 5(j) in the large aspect-ratio approximation (52) and such 
an arbitrary phase is still present in the unapproximated leading order DC equations, as 
discussed in Sec. [| However, the full numerics show that there is a slow drift towards a 
specific set of 8(j) for each pattern. Several authors |7], ||, |10[ have studied this inter-row 
phase locking in Pattern 1, but a satisfying answer is yet to be developed. The zero-field 
limit (/ = 0) is an exception in that exact neutral stability and a family of periodic solutions 
can be found || |33|, implying that there is no inter- row locking. On the other hand, a slow 



drift starts to occur as / is perturbed away from zero fl4"4" |. We conjecture that the arbitrary 
phases should be constrained by a solvability condition in the higher-order expansions of our 
analysis, which is automatically satisfied when / = 0. Finding that condition, however, is 
likely to be an elaborate task. 

Our analysis is based on such simplifications as zero temperature, no disorder, and no 
self-fields. Clearly, the effect of relaxing these assumptions should be also investigated. 
Thermal noise, self-fields and inhomogeneities alter the switching sequence in simulations 
of the row-switching cascade |22j, [32], |18j. This might explain the irregularity of the row- 
switching order observed experimentally by Trias and Lachenmann et al. |20|. On the 
other hand, the directed use of disorder (e.g., by removing some of the edges in the array) 
might prove a valuable strategy to enhance the locking property of the arrays [49|. Including 
inductances would also change the current distributions j|, [31], [32]]. Previous work [^, |TE| , 
and our own preliminary calculations including self-inductances, show that RS states persist 
at least for small inductances. Our expansion could be extended to include inductances and 
then proceed to describe the modified solutions. However, qualitatively new phenomena can 
also arise. For example, it is known |19, [26|] that 



if any inductance is included in the 
model, a coherent state (dynamical checkerboard pattern) emerges near / = 1/2 when the 
RS states cease to exist. 

In this article, we have only considered "clean" RS states, formed by whirling and no- 
vortex superconducting regions. Thus, we have assumed that the Q rows do not contain any 
static vortices. It is generally expected that the depinning of a Q region would become easier 
when it contains a pre-trapped vortex. Therefore, the existence of states with static vortices 
probably does not affect the critical curve in Fig. H](c). However, the question of how the 
depinning of a static 2D array depends on various parameters (T, /, N x and N y ) is not fully 
understood, except in the case of the ladder [J7|, and requires further scrutiny. Similarly, 
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Figure 15: 
size N x = 
Fig 



11 



Dimensionless AC voltage amplitudes <p y for Pattern 1 in a large array of 
= 63 and N y = 31. The other parameters are the ones used in Pattern 1 in 
The amplitudes decay quickly from the boundaries, and nearly vanish inside the 
array. The tendency was already present in Fig. [TTJ for the 31 x 7 array, but it is clearer 
here. Consequently, were it to be used as an oscillator, the total AC output voltage for this 
pattern would not scale favorably with increasing size. 



the S rows in the RS states were assumed to be in the whirling (normal resistive) state. Our 
simulations sometimes show "generalized" RS states which contain one or more rows that are 
neither switched nor quiescent, but "active". The states could be born, for instance, when 
Idc is increased so that vortices start to enter a Q region but not strongly enough to switch it. 
Junctions in the active rows undergo highly nonlinear oscillations, and propagating vortices 
are localized. These states create additional steps in the I—V characteristics between two RS 
steps, and are detectable. Thus, they should be considered for a comprehensive treatment 
of row-switching. 

Apart from investigating the RS states, we have introduced in this article a systematic 
approach to the analysis of the dynamics of 2D Josephson arrays. Unlike ID arrays, which 
have already led to a great amount of insight into important phenomena (such as soliton 
propagation and interaction in the parallel-connected arrays |5(J, or synchronization, 
clustering, and magnet-like phase transitions in the series-connected arrays [0, fjO|), 2D 
arrays have been much harder to analyze. This is partly due to their network equations 
being more complicated, and also to their having a wider variety of solutions. As our weakly- 
nonlinear analysis shows, the difficulty regarding the formulation is reduced by the compact 
mesh formalism introduced in the previous numerical studies [^7], ^9], f|, [31], |32|, |18R . We 
feel that the transparent form of the mesh equations has the potential to provide analytical 
information in the strongly nonlinear regime. 

Of these strongly nonlinear solutions, two are of particular interest. First, coherent states, 
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such as Pattern 1 in Fig. ^], might be suitable for oscillator applications, if / is kept small 
so that the whole array operates nearly in-phase. However, for the completely row-switched 
state to be useful, overdamped junctions, which rotate less smoothly and, thus, produce 
larger AC amplitudes, should be employed ||. The extension of our analysis for this case 
(concerning only Pattern 1) appears to be straightforward. However, we can already point 
out a complication due to the spatial distribution of the AC amplitudes. Recall how the 
<py amplitudes in Pattern 1 (Fig. |ll|) decay from the boundaries into interior of the array. 
This effect is more clearly illustrated in Fig. 15 where the AC amplitudes are computed for 
Pattern 1 in a larger array (N x = 63 and N y = 31), the other parameters being identical. 
The amplitudes decay quickly, and nearly vanish inside the array. Consequently, the total 
AC voltage does not increase significantly even when more junctions are inter-connected. 

Finally, flux flow [37], |32], |I^| is also a highly nonlinear but disordered regime in which 
localized vortices propagate "diffusively". Theoretical studies so far have been based on 
phenomenological pictures of vortices and their interactions . A more formal treatment of 
these solutions and a detailed prediction of, for instance, the flux flow resistance is strongly 
awaited both from the theoretical and experimental points of view. 
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A Small aspect-ratio approximation 



Following the large aspect-ratio approximation presented in Section |3.2.1| , we calculate now 
a bulk approximation to the DC equations in a Q region with a small aspect- 

ratio, i.e. when its vertical size N y is much larger than N x . Far from the top and bottom 
boundaries, the solution is expected to be independent of j (assuming there are no trapped 
vortices). Then, the DC equations (|29|),(|30|) simplify to: 



sin 



o 



thus 



o 



(65) 



W -r >-j ~ vW) = ( 66 ) 

From ( |52"D and the boundary conditions fl3"7] ) we can construct the following telescope sum 
which must be satisfied 



£ 

i=i 



sin ■ 



(N x + 1) J dc . 



(67) 



From these two equations (EH) , (|67D we can then solve for the vertical phases in the bulk of 
the Q region: 

0K(t) = 2tt/ + l-ij+a (68) 

where 



a 



. [ (N* + 1) smjnf) ' 
arcsm < r , „ — —Idc 



sm[(N x + 1)tt/] 
Compare this with the large aspect-ratio case fl47|) in which 



= arcsin Jd c is independent of 
% in the bulk of a Q region. In contrast, in the present small aspect-ratio case, the external 
field / is absorbed now by the vertical junctions in order to ensure the flux quantization 
restriction Note also how, in this case, consideration of the top and bottom edges, 

neglected from the bulk Q region, becomes crucial to introduce matching across the switched 
regions or to the array boundaries. Without the correction from the edges, phase relations 
across the S rows are not well defined. However, the small aspect-ratio approximation is still 
significant because it provides a clue to an important question: what is the lower bound for 
a Q region to remain unbroken? Thus, we use this calculation when we discuss the existence 
and stability of RS patterns in Sec. || In this context, the small aspect-ratio approximation 
is the limiting case for which a Q region is most easily broken by raising either / or I^ c . 
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